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LONG-TERM  GOALS 

The  long-term  goal  is  to  develop  a  nonhydrostatic,  parallel  ocean  simulation  tool  that  is  capable  of 
simulating  processes  on  a  wide  range  of  scales  through  use  of  accurate  numerical  methods  and  high- 
performance  computational  algorithms.  The  tool  will  be  applied  to  study  highly  nonlinear  internal 
waves  in  coastal  domains  to  develop  an  improved  understanding  of  mechanisms  that  govern  their 
generation,  propagation,  and  dissipation. 

OBJECTIVES 

The  primary  objective  is  to  enhance  the  capabilities  of  the  SUNTANS  model  through  development  of 
algorithms  to  study  multiscale  processes  in  estuaries  and  the  coastal  ocean.  This  involves  development 
of  1)  improved  momentum  and  scalar  advection  on  unstructured,  staggered  grids,  2)  accurate  and 
efficient  algorithms  for  solution  of  the  nonhydrostatic  pressure,  and  3)  adaptive  grid  capabiliites  with 
adaptive  mesh  refinement  and  model  nesting. 

APPROACH 

This  work  focuses  on  the  continued  development  of  SUNTANS  (Stanford  Unstructured 
Nonhydrostatic  Terrain-following  Adaptive  Navier-Stokes  Simulator),  a  free-surface,  nonhydrostatic, 
unstructured-grid,  parallel  coastal  ocean  and  estuary  simulation  tool  that  solves  the  Navier-Stokes 
equations  under  the  Boussinesq  approximation  (Fringer  et  al.,2006).  The  formulation  is  based  on  the 
method  outlined  by  Casulli  and  Walters  (2000),  in  which  the  free-surface  and  vertical  diffusion  are 
discretized  with  the  theta  method  which  eliminates  the  Courant  condition  associated  with  fast  free- 
surface  waves  and  the  elevated  friction  term  associated  with  small  vertical  grid  spacing  at  the  free- 
surface  and  bottom  boundary.  For  flows  with  extensive  wetting  and  drying,  advection  of  momentum  is 
accomplished  with  the  semi-Lagrangian  advection  scheme  (Wang  et  al.  201  la),  which  ensures  stability 
in  the  presence  of  cells  that  fill  and  empty  with  the  tides.  Scalar  advection  is  accomplished  semi- 
implicitly  and  continuity  of  volume  and  mass  are  guaranteed  for  the  hydrostatic  solver.  The  theta 
method  for  the  free  surface  yields  a  two-dimensional  Poisson  equation,  and  the  nonhydrostatic 
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pressure  is  governed  by  a  three-dimensional  Poisson  equation.  These  are  both  solved  with  the 
preconditioned  conjugate  gradient  algorithm  with  Jacobi  and  block-Jacobi  preconditioning, 
respectively.  Because  the  nonhydrostatic  component  of  SUNTANS  is  essentially  a  correction  to  the 
hydrostatic  component,  SUNTANS  can  be  run  seamlessly  in  nonhydrostatic  or  hydrostatic  modes. 
SUNTANS  is  written  in  the  C  programming  language,  and  the  message-passing  interface  (MPI)  is 
employed  for  use  in  a  distributed-memory  parallel  computing  environment.  SUNTANS  employs  the 
generalized  length  scale  approach  to  Reynolds-averaged  turbulence  modeling  (Wang  et  al.  2011b). 

The  SUNTANS  grid  employs  z-levels  in  the  vertical  and  is  unstructured  in  plan,  which  enables  the 
resolution  of  complex  coastlines  and  topographic  features.  Unstructured  grids  also  enable  the  use  of 
high  grid  resolution  in  regions  of  interest  while  coarsening  the  grid  in  regions  where  grid  resolution  is 
not  required,  thereby  significantly  reducing  computational  overhead. 

WORK  COMPLETED 

We  completed  the  development  of  a  nonhydrostatic  isopycnal-coordinate  model  (Vitousek  and  Fringer 
2014)  and  here  present  results  of  the  test  cases  from  that  work  (that  were  not  reported  in  last  year’s 
annual  report).  Follow-on  work  as  part  of  the  FLEAT  and  NOPP  DRIs  will  focus  on  implementing  the 
isopycnal-coordinate  model  into  the  SUNTANS  modeling  framework.  We  have  also  performed  large- 
eddy  simulations  of  breaking  internal  waves  on  slopes  and  computed  the  associated  cross-shelf 
transport.  Although  these  simulations  do  not  employ  the  SUNTANS  model,  it  is  our  hope  that  the 
SUNTANS  model  will  eventually  have  this  capability. 

RESULTS 

Nonhydrostatic  isopycnal-coordinate  ocean  model 

We  completed  a  manuscript  reporting  the  formulation  and  testing  of  a  nonhydrostatic  isopycnal 
coordinate  ocean  model  (Vitousek  and  Fringer  2014),  which  to  our  knowledge  represents  the  first  of  its 
kind.  While  numerous  ocean  models  are  written  in  isopycnal  coordinates  (most  notably 
MICOM/HYCOM;  Bleck  et  al.,  1992;  Bleck,  2002),  our  model  is  unique  because  it  solves  for  the 
nonhydrostatic  pressure.  This  enables  us  to  compute  nonhydrostatic  internal  solitary  waves  at  a 
fraction  of  the  computational  cost  of  traditional  z-coordinate  models  (e.g.  Zhang  et  al.  2011)  because 
the  number  of  vertical  coordinates  can  be  reduced  from  0(100)  to  0(1)  while  still  resolving  a  bulk  of 
the  internal  wave  energetics.  Below  we  report  on  two  test  cases  that  highlight  the  features  of  the 
model.  More  detail  on  these  test  cases  can  be  found  in  the  paper  by  Vitousek  and  Fringer  (2014). 

Formation  of  internal  solitary  waves  in  a  laboratory  experiment:  This  test  case  compares  the  formation 
of  an  internal  solitary-like  wave  train  using  the  isopycnal-coordinate  model  to  the  laboratory 
experiments  in  Horn  et  al.  (2001).  The  laboratory  experiments  investigate  the  formation  of  solitary 
waves  in  an  approximately  two-layer  density  profile  with  a  pycnocline  thickness  of  0.01  m  in  an 
enclosed  tank  of  depth  0.29  m  and  length  6  m.  The  internal  interface  is  initialized  with  a  linear  tilt 
with  a  wave  amplitude  given  by  [0.01305,  0.0261,0.03915,  0.0522,  0.0783]  m.  The  model  of  the 
laboratory  experiments  is  run  with  2  layers  and  256  horizontal  grid  points  and  with  a  slightly  lower 
density  difference  than  the  lab  experiment  because  the  finite-thickness  pycnocline  in  the  laboratory 
experiment  reduces  the  wave  speed  slightly  when  compared  to  the  two-layer  model  configuration.  A 
drag  coefficient  Cd=0.0175  is  employed,  and  the  horizontal  and  vertical  eddy-viscosities  are  equal  and 
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given  by  10'  m  s'  .  These  parameters  are  chosen  to  fit  the  amplitude  of  the  solitary  waves  formed 


during  the  laboratory  experiment.  Figure  1  shows  the  results  of  the  model-predicted  time  series  of  the 
interface  displacements  (located  at  the  center  of  the  tank)  compared  to  the  laboratory  experiments  of 
Horn  et  al.  (2001).  Agreement  between  the  model  and  the  laboratory  experiments  is  excellent  for  the 
first  half  of  the  experiment  (t  <  200  s).  However,  for  the  second  half  of  the  experiments  there  is  some 


mismatch  in  the  phase  and  amplitude  of  the  waves.  Over  the  course  of  a  laboratory  experiment  Horn  et 


al.  (2001)  reported  “a  gradual  thickening  of  the  density  interface  over  the  set,  typically  from 
approximately  1  cm-2  cm”.  The  finite-thickness  interface  and  thickening  are  partly  responsible  for  the 
mismatch  with  the  model  (which  has  a  sharp  interface  and  no  interface  thickening).  Additionally,  the 
model  appears  to  form  slightly  wider  solitary  waves  than  the  laboratory  experiment.  This  indicates  that 
the  dispersive  properties  of  the  waves  are  not  captured  exactly  in  the  model  because  of  the  two-layer 
approximation.  Overall,  however,  the  model  adequately  captures  the  formation  of  solitary  waves  in  the 
laboratory  experiment. 


- Nonhydrostatic  isopycnal  model 

- Horn  2001 


5 


[a]  a  =  1.305  cm 


-5 

5 


\b\  a.  =  2.61  cm 


-5 


5m — T - 

[c]  a.  =  3.915  cm 


-5 


\~D\  a  =  5  7?  cm 


-5 


-5 


0 


50  100  150  200  250  300  350  400 

time  [s] 


Figure  1:  Comparison  of  nonhydrostatic  model-predicted  and  measured  interface  displacement  time 
series  at  the  middle  of  the  tank  from  the  experiments  of  Horn  et  al.  (2001). 

Nonhydrostatic  internal  wave  generation  over  a  Gaussian  ridge:  Here  we  simulate  the  generation  of 
internal  tides  over  Gaussian  ridge  topography  following  the  paper  by  Buijsman  et  al.  (2010).  The 
model  is  run  with  bathymetry  that  approximates  the  eastern  ridge  of  the  Luzon  Strait  with  a  height  of 
2600  m  in  a  depth  of  3000  m,  as  shown  in  Figure  2.  This  domain  is  800  km  long,  and  1750  grid  points 
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are  used  in  the  horizontal,  which  ensures  that  leading-order  nonhydrostatic  effects  are  resolved 
(Vitousek  and  Fringer  2011).  The  boundaries  are  forced  with  barotropic  currents  oscillating  at  the  M2 
tidal  frequency  and  with  a  magnitude  of  0.134  m  s'1.  The  initial  2-  and  10-layer  configurations  are  set 
to  best  match  the  weakly  nonlinear  and  nonhydrostatic  behavior  of  solitary  waves  in  the  stratification 
of  Buijsman  et  al.  (2010),  which  is  representative  of  the  South  China  Sea.  The  still-water  isopycnal 
depths  in  the  10-layer  model  are  set  to  z=0,  50,  100,  250,  460,  750,  1000,  1500,  2000,  2500,  3000  m 
with  density  values  computed  at  the  vertical  center  of  each  layer.  The  still-water  internal  interface  for 
the  2-layer  model  is  set  to  460  m  and  the  upper-  and  lower-layer  densities  are  set  to  match  the  mode-1 
wave  speed  of  the  10-layer  model.  A  horizontal  eddy- viscosity  of  100  m2  s'1  and  a  vertical  eddy- 
viscosity  of  10"  m  s'  are  used  to  characterize  unresolved  mixing  of  momentum  due  to  turbulence  in 
the  presence  of  internal  waves,  along  with  a  bottom  drag  coefficient  of  Cd=0.0025. 

Figure  3  compares  the  generation  and  evolution  of  internal  solitary  waves  using  2  and  10  layers. 
Although  both  models  form  trains  of  nonlinear  internal  solitary-like  waves  with  nearly  identical  wave 
speeds,  the  amplitude,  length  scales,  and  number  of  solitary  waves  differ  between  the  models.  The 
two-layer  model  only  supports  first-mode  waves,  whereas  the  10-layer  model  supports  higher  modes. 
Consequently  barotropic-to-baroclinic  conversion  in  the  2-layer  model  is  restricted  to  the  first 
baroclinic  mode,  whereas  energy  is  converted  to  higher-mode  waves  in  the  10-layer  model.  This 
results  in  larger  amplitude  solitary  waves  in  the  2-layer  model  relative  to  the  10-layer  model. 
Additionally,  the  solitary  wave  widths  in  the  2-layer  model  are  slightly  larger  than  those  in  the  10-layer 
model  due  to  unresolved  nonhydrostatic  dispersion  in  the  2-layer  model. 
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Figure  2:  Model  domain  and  bathymetry  used  in  the  South  China  Sea  case  following  Buijsman  et 
al.  (2010).  The  black  filled  area  represents  the  bathymetry,  the  blue  line  represents  the  free  surface, 
and  the  black  dashed  lines  represent  the  locations  of  the  sponge  layers. 
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Figure  3:  Comparison  of  internal  wave  generation  in  the  nonhydrostatic  2-layer  (red)  and  10-layer 
(blue)  models  at  time  intervals  indicated  by  the  red  dot  in  the  barotropic  current  time  series  in  the 

panels  in  the  right  column. 


Transport  by  breaking  internal  waves  on  slopes 

Using  the  code  of  Cui  (1999),  Arthur  and  Fringer  (2014)  performed  three-dimensional  DNS  of 
breaking  solitary-like  internal  waves  on  slopes  to  examine  the  energetics  and  mixing  as  a  function  of 
internal  wave  steepness  and  bathymetric  steepness.  Here  we  report  on  the  results  of  transport  by 
breaking  internal  waves  computed  with  this  same  setup.  Figure  4  shows  snapshots  in  time  of  a  plume 
of  particles  initially  located  where  the  pycnocline  intersects  the  slope.  As  the  wave  begins  to  interact 
with  the  slope,  particles  are  transported  offshore  in  the  lower  layer  and  slightly  onshore  in  the  upper 
layer  (Figure  4a-c).  After  the  wave  breaks,  the  upslope  surge  of  dense  fluid  transports  particles 
onshore.  Upper-layer  fluid  is  directed  up  and  over  the  nose  of  the  upslope  surge,  carrying  a  thin  layer 
of  particles  offshore  along  the  pycnocline  (Figure  4d-e).  Once  the  surge  reaches  its  maximum  upslope 
location,  the  relaxation  of  dense  fluid  transports  particles  back  downslope  (Figure  4f).  When  the  jet 
detaches  from  the  slope  as  an  intrusion,  particles  are  transported  further  offshore  along  the  pycnocline 
(Figure  4f-i). 
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Over  the  course  of  the  breaking  event,  onshore  transport  occurs  along  the  bottom  within  the  upslope 
surge  of  dense  fluid.  Offshore  transport  occurs  along  the  pycnocline,  initially  due  to  flow  over  the  nose 
of  the  upslope  surge,  and  ultimately  due  to  the  offshore  intrusion  of  intermediate  density  fluid.  Above 
the  pycnocline,  there  is  minimal  onshore  transport  that  increases  with  depth.  The  shape  of  the  plume  at 
the  end  of  the  breaking  event  resembles  the  plume  in  the  numerical  dye  study  of  Nakayama  and 
Imberger  (2010).  Their  study  examined  the  cumulative  effect  of  multiple  periodic  internal  waves 
breaking  on  a  slope,  as  opposed  to  the  single  wave  of  depression  considered  here. 
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Figure  4:  Snapshots  in  time  of  a  plume  of particles  for  a  breaking  internal  wave.  The  cross-shore 
position  is  presented  relative  to  the  initial  mean  position  of  the  plume  xf°and  normalized  by  the 
solitary  wave  length  Lir  The  pycnocline  is  shown  in  black,  while  the  top  panel  (labeled  “ Initial ”) 

shows  the  initial  particle  plume. 

The  layer  of  particles  transported  offshore  along  the  pycnocline  by  the  breaking  wave  in  Figure  4 
resembles  an  intermediate  nepheloid  layer  (INL;  McPhee-Shaw  et  al.  2004,  McPhee-Shaw  2006, 
Cheriton  et  al.  2014).  If  we  define  the  INL  as  the  group  of  particles  within  the  density  isosurfaces 
comprising  the  pycnocline,  then  the  simulations  can  be  used  to  determine  the  initial  location  of  these 
particles.  Figure  5  shows  the  INL  formed  at  the  end  of  the  breaking  event  depicted  in  Figure  4,  along 
with  the  location  of  the  particles  that  comprise  the  INL  prior  to  the  breaking  event.  This  figure  shows 
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that  mixing  by  the  breaking  internal  wave  causes  particles  to  be  entrained  from  above  and  below  the 
pycnocline  into  the  INL.  In  this  particular  case,  92%  of  the  particles  that  comprise  the  INL  originate 
outside  of  the  pycnocline. 


Figure  5:  Particles  that  make  up  the  intermediate  nepheloid  layer  for  the  breaking  wave  shown  in 
Figure  4.  (a)  Shows  the  final  positions  of  the  particles  within  the  INL  at  the  end  of  the  breaking 
event  (time  t  =  tf  while  (b)  shows  their  initial  positions  (time  t  =  t0).  Each  plot  shows  density 
isosurfaces  that  make  up  the  pycnocline  as  well  as  a  vertical  line  at  the  cross-shore  position  xltPyc  of 

the  intersection  of  the  initial  pycnocline  and  the  slope. 

The  lateral  turbulent  velocities  induced  by  breaking  in  the  three-dimensional  simulations  lead  to  lateral 
transport  of  particles,  as  shown  in  Figure  6.  Before  turbulence  develops,  lateral  transport  occurs  only 
due  to  a  random  walk  based  on  the  molecular  diffusivity.  However,  after  breaking,  the  flow  becomes 
turbulent  and  lateral  transport  occurs  at  a  faster  rate.  This  increased  lateral  spreading  is  evident  during 
the  upslope  surge  (Figure  6e),  and  continues  for  the  rest  of  the  breaking  event  (Figure  6f-g).  Most 
lateral  transport  appears  to  occur  during  the  downslope  flow  of  dense  fluid  after  the  upslope  surge  and 
the  resulting  intrusion  of  intermediate  density  fluid  along  the  interface.  Note  that  in  Figure  6,  particles 
are  plotted  outside  of  the  lateral  boundaries  of  the  computational  domain.  This  is  because  lateral 
particle  positions  were  adjusted  to  account  for  periodic  boundary  crossings  by  adding  or  subtracting 
the  width  of  the  domain. 
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Figure  6:  Lateral  transport  of  particles  for  a  breaking  internal  wave.  Each  panel  shows  a  top  view 
(xrx2)  of  all  particles.  For  clarity,  particles  are  shown  to  extend  beyond  the  lateral  boundaries  of 
the  domain  owing  to  lateral  periodicity.  Also  included  are  lines  representing  three  standard 
deviations  of  the  lateral  particle  positions,  which  indicate  the  lateral  boundaries  of  the  particle 
plume.  Labels  d-g  correspond  to  the  time  snapshots  shown  in  Figure  4. 


IMPACT/APPLICATIONS 


High-resolution  simulations  using  nonhydrostatic  models  like  SUNTANS  are  crucial  for  understanding 
multiscale  processes  that  are  unresolved,  and  hence  parameterized,  in  larger-scale  ocean  models. 
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